%   brk2msk.m - outputs the mri scalp mask as both a volume mask and a list of coord
%   Chris Whalen
%   University of Illinois
%   06/24/04 - Date Created
%   08/19/05 - Last Modified


function [xyz_mask, vol_mask, vol] = brk2msk(file,mri_dim,vol,thresh)

% fprintf('\n\n     brk2msk \n');

% tic

% Reading and Initializing vol=0
if vol == 0
    fid = fopen(file.brik);
    data=fread(fid,'uint16');
    fclose(fid);
    data=data/(max(data));  % rescale intensities 0 to 1
    vol = reshape(data,[mri_dim.x mri_dim.y mri_dim.z]);    % create BRIK
    vol = filter3d(vol); % Spatial Filters to remove graininess 
    vol = smooth3(vol); % gaussian filtering with kernel 3 and sd 0.65 . Must be done here or else xyz_mask won't be smoothed
    clear data
end

% Initializing
vol_mask = zeros(mri_dim.x, mri_dim.y, mri_dim.z); 
max_msk_pts = 192000;             % Maximum Points that could be contained in Mask
xyz_mask = zeros(max_msk_pts,3);  % This is over the maximum number of entries possible; Preallocating speeds up processing
temp_mask = zeros(max_msk_pts,3); % Preallocating speeds up processing

% % Thresholding
% fprintf('Enter Scalp Threshold (0-1):  (Default Value 0.05)\n');
% thresh = input ('Threshold: ', 's');
% if strcmp(thresh,'') ==1
%   thresh = 0.05;
% else
%   thresh = str2num(thresh);
% end

% File Being Processed
% fprintf('\nSegmenting Scalp for %s ... \n', file.brik);
% fprintf('\nEstimated Segmentation Time: 0.2 minutes\n');
% tic;
  
% Generating Mask Left/Right
% fprintf('Generating L-R Mask ...  \n');
%fprintf('Generating L-R Mask ... \n  Estimated Time: 0.2 min \n');
%tic;
change = 0;
mask_row = 1;                                   
for k=1:mri_dim.z                                       % Z Coords
    for j=1:mri_dim.y                                   % Y Coords
        for i=1:mri_dim.x                         % X Right Bound Scalp Coords (must do at least slightly more than half way in case head is asymetrical)
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 xyz_mask(mask_row,:) = [i j k];
                 mask_row = mask_row + 1;
                 change = 1;
                 vol_mask(i,j,k) = 1;
            end               
        end
        change = 0;
        for i=mri_dim.x:-1:(1)             % X Left Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 xyz_mask(mask_row,:) = [i j k];
                 mask_row = mask_row + 1;
                 change = 1;
                 vol_mask(i,j,k) = 1;
            end               
        end
        change = 0;    
    end
end
%t_left_right = toc;
%fprintf('  Actual Time: %-2.1f min\n', t_left_right/60);

% Generating Mask Anterior/Posterior
% fprintf('Generating A-P Mask ... \n');
%fprintf('Generating A-P Mask ... \n  Estimated Time: 0.2 min \n');
%tic;
change = 0;                                 
for k=1:mri_dim.z                                       % Z Coords
    for i=1:mri_dim.x                                   % Y Coords
        for j=1:(mri_dim.y)                         % Anterior Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 xyz_mask(mask_row,:) = [i j k];
                 mask_row = mask_row + 1;
                 change = 1;
                 vol_mask(i,j,k) = 1;
            end               
        end
        change = 0;
        for j=(mri_dim.y):-1:1             % Posterior Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 xyz_mask(mask_row,:) = [i j k];
                 mask_row = mask_row + 1;
                 change = 1;
                 vol_mask(i,j,k) = 1;
            end               
        end
        change = 0;    
    end
end
%t_ant_post = toc;
%fprintf('  Actual Time: %-2.1f min\n', t_ant_post/60);

% Generating Mask Superior/Inferior
% fprintf('Generating S-I Mask ... \n');
%fprintf('Generating S-I Mask ... \n  Estimated Time: 0.2 min \n');
%tic;
change = 0;                           
for i=1:mri_dim.x                                       % Z Coords
    for j=1:mri_dim.y                                   % Y Coords
        for k=(mri_dim.z):-1:1             % Superior Bound Scalp Coords 
            if  (vol(i,j,k) >= thresh) & (change == 0) %(sum(vol(i,j,k)) == lin_clust_size)
                 xyz_mask(mask_row,:) = [i j k];
                 mask_row = mask_row + 1;
                 change = 1;
                 vol_mask(i,j,k) = 1;        
            end               
        end
        change = 0;    
    end
end
%t_sup_inf = toc;
%fprintf('  Actual Time: %-2.1f min\n', t_sup_inf/60);

xyz_mask(mask_row:max_msk_pts,:) = []; % Reducing extra entries

% Sort and Remove Identical Scalp Voxels
% fprintf('Condensing Mask Data Set... \n');
%fprintf('  Estimated Time: 0-1 min\n');
xyz_mask = sortrows(xyz_mask,[1 2 3]);
temp_row = 1;
%tic;
for row=1:((size(xyz_mask, 1))-1)
  if (xyz_mask(row,1) == xyz_mask(row+1,1)) & (xyz_mask(row,2) == xyz_mask(row+1,2)) & (xyz_mask(row,3) == xyz_mask(row+1,3))  
  else 
      temp_mask(temp_row,:) = xyz_mask(row, :);
      temp_row = temp_row + 1;
  end
end
temp_mask(temp_row,:) = xyz_mask((size(xyz_mask, 1)),:);
temp_mask((temp_row+1):max_msk_pts,:) = []; % Delete excess preallocated rows
rep_pts = size(xyz_mask,1) - size(temp_mask,1);
xyz_mask = temp_mask;
%t_reduce_mask = toc;

%fprintf('  Actual Time: %-3.0f min\n', t_reduce_mask/60);
% fprintf('Number of repeated points: %-3.0f\n', rep_pts);
% fprintf('Number of  scalp   points: %-3.0f\n', length(xyz_mask));

% xyz_mask = xyz_mask + 1; % Adding 1 to have both BRIK and ELP be in same dimensions
% t_brk2msk = toc;
% fprintf('  Actual Time: %-3.0f min\n', t_brk2msk/60);
% 
% fprintf('\n     brk2msk Done \n');





